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Abstract 

We report the coupled effects of macroion charge discretization and counterion va- 
lence in the primitive model for spherical colloids. Instead of considering a uniformly 
charged surface, as it is traditionally done, we consider a more realistic situation 
where discrete monovalent microscopic charges are randomly distributed over the 
sphere. Monovalent or multivalent counterions ensure global electroneutrality. We 
use molecular dynamics simulations to study these effects at the ground state and 
for finite temperature. The ground state analysis concerns the counterion structure 
and charge inversion. Results are discussed in terms of simple analytical models. 
For finite temperature, strong and weak Coulomb couplings are treated. In this sit- 
uation of finite temperature, we considered and discussed the phenomena of ionic 
pairing (pinning) and unpairing (unpinning). 
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1 Introduction 



Charged colloidal suspensions are a subject of intense experimental and the- 
oretical work not only because of their direct application in industrial or bio- 
logical processes, but also because they represent model systems for atomistic 
systems. The electrostatic interactions involved in such systems have a funda- 
mental role in determining their physico-chemical properties [1,2]. Theoretical 
description of highly charged colloidal solutions faces two challenges: (i) differ- 
ent typical length scales due to the presence of macroions (i. e. charged colloids 
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of the size 10-1000 A) and microscopic counterions and (ii) their long-range 
Coulomb interaction. A first simphfying assumption is to treat the solvent as 
a dielectric medium solely characterized by its relative permittivity e^. A sec- 
ond widely used approximation consists in modeling the short range ion-ion 
excluded volume interaction by hard spheres. These two approximations are 
the basis of the so-called primitive model of electrolyte solutions. The system 
under consideration is an asymmetrical electrolyte solution made up of highly 
charged macroions and small counterions. A further simplification motivated 
by this asymmetry can be made by partitioning the system into subvolumes 
(cells) , each containing one macroion together with its neutralizing counterions 
plus (if present) additional salt. This approximation has been called accord- 
ingly the cell model [3,4]. The cells adopt the symmetry of the macroion, here 
spherical, and are electrostatically decoupled. It is within the cell model that 
we present our simulation results. 

For spherical macroions the structural charge is usually modeled by a central 
charge, which by Gauss theorem is equivalent to a uniform surface charge 
density as far as the electric field (or potential) outside the spherical colloid 
is concerned. 

Most analytical concepts as well as simulations rely on the above assumptions 
and especially on the central charge assumption. It is well known that in the 
strong Coulomb coupling regime ion-ion correlations become very important, 
and significant deviations from mean-field approaches are expected. A counter- 
intuitive effect which classical mean-field theories (like Poisson-Boltzmann 
model) cannot explain is the phenomenon of overcharge, also called charge in- 
version. That is, there are counterions in excess in the vicinity of the macroion 
surface so that its net charge changes sign. This has recently attracted signifi- 
cant attention [5-17]. In particular, we showed recently that this phenomenon 
may give rise to a strong long range attraction between like-sign charged col- 
loids [12,13,17]. A natural question which comes up is: does overcharge and 
more generally ion-ion correlations strongly depend on the way the macroion 
structural charge is represented (i. e. uniformly charged or discrete charges on 
its surface)? In a recent paper [16], we studied such a situation in the strong 
Coulomb coupling regime where the macroion charge was carried by divalent 
microions in the presence of divalent counterions {same ionic valence). In Ref. 
[16] we reported the important result showing that overcharge is still pos- 
sible under those conditions. Moreover we showed that the intrinsic electric 
field solely due to the macroion surface microions (without counterions) varies 
strongly from point to point on the colloidal sphere [16]. 

The goal of this paper is to study by means of molecular dynamics (MD) sim- 
ulations the coupled effects of macroion charge discretization and counterion 
valence in the primitive model for spherical colloids. A systematic comparison 
with the uniform macroion charge distribution (i. e. central charge) is under- 
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taken. The paper is organized as follows. In Sec. 2 we give some details on the 
macroion charge discretization as well as on the MD simulation model. Section 
3 is devoted to the ground state analysis where surface counterion structure 
and overcharge are addressed. In Sec. 4 we investigate the finite temperature 
situation, where counterion structure is studied for strong and weak Coulomb 
couplings. Finally, in Sec. 5 we provide a summary of the results. 



2 Simulation model 

2.1 Macroion charge discretization 

The procedure is similar to the one used in a previous study [16]. The discrete 
macroion charge distribution is produced by using Zm monovalent microions 
of diameter a (same diameter as the counterions) distributed randomly on the 
surface of the macroion. Then the structural charge is Q = —ZmQ — —Z^ZdC, 
where Zm > 0, Z^, — 1 is the valence of these discrete microions and e is the 
positive elementary charge. These discrete colloidal charges (DCC) are fixed 
on the surface of the spherical macroion. Figure 1 shows a schematic view of 
the setup. The counterions (not shown in Fig. 1) have a charge q = +Zce, 
where Z^ > stands for the counterion valence. In spherical coordinates the 
elementary surface is given by: 

dA — Tq sin 9d9dip = — roC?(cos 9)dip , (1) 

where Tq is the distance between the macroion center and the DCC center. 
Thus to produce a random discrete charge distribution on the surface we gener- 
ated (uniformly) randomly the variables cos 9 and Lp. Excluded volume is taken 
into account by rejecting configurations leading to an overlap of microions. 
Phenomena such as surface chemical reactions [18], hydration, roughness [19] 
are not considered. For commodity we introduce the notation {—Z^ '■ +Zc) to 
characterize the valence symmetry (asymmetry) for Zc = 1 {Zc > 1) of the 
ions (DCC and counterions) involved in discrete systems. 



2.2 Molecular dynamics procedure 

A MD simulation technique was used to compute the motion of the counterions 
coupled to a heat bath acting through a weak stochastic force W(t). The 
procedure is very similar to the one used in previous studies [12,16]. 
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Fig. 1. Schematic view of the setup: the discrete colloidal charges (DCC) of diameter 
a are in dark grey. For a detailed meaning of the other symbols see text. Note that 
this a a two-dimensional representation of the three-dimensional system. 

The motion of counterion i (DCC ions being fixed) obeys the Langevin equa- 
tion 

= -V.t/(r.) - + W,(t) , (2) 



where m is the counterion mass, U is the potential force having two contri- 
butions: (i) the Coulomb interaction and (ii) the excluded volume interaction, 
and 7 is the friction coefficient. Friction and stochastic force are linked by the 
dissipation- fluctuation theorem < Wj(t) ■ Wjit') >= 6m'ykBT6ij6(t — t). For 
the ground state simulations the stochastic force is set to zero. 

Excluded volume interactions are taken into account with a pure repulsive 
Lennard- Jones potential given by 



ULj(r) 




( — ) - ( — ) 



for r — ro < r, 



cut: 



for r - ro > r^ut, 



(3) 



where Tq = for the microion-microion interaction (the microion being either 
a counterion or a DCC), ro = 7a for the macroion-counterion interaction 
and Tcut = 2^/^(7 is the cutoff radius. This leads to a (center-center) macroion- 
counterion distance of closest approach a = 8a (see also Fig. 1). The macroion 
surface charge density am is defined as 
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Table 1 

Simulation parameters with some fixed values. 



parameters 


a = 3.57 A 


Lennard Jones length units 


To = 298K 


room temperature 


eij = ksTo 


Lennard Jones energy units 


Zjn 


macroion valence 




discrete colloidal charge valence 




counterion valence 


Ib 


Bjerrum length 


a = 8a 


macroion-counterion distance of closest approach 




macroion surface charge density 


R = 40a 


simulation cell radius 


fm = 8x 10-3 


macroion volume fraction 



Energy and length units in our simulations are related to experimental units 
by taking eij =A;bTo (with Tq = 298 K) and a = 3.57A respectively. 

The pair electrostatic interaction of any pair ij, where i and j denote either a 
DCC a counterion or the central charge (for the non-discrete case), reads 

Ucouiir) = kBToh^ , (5) 



where Is = e"^ /ATreoerksTo is the Bjerrum length describing the electrostatic 
strength. For the rest of this paper, electrostatic energy will always be ex- 
pressed in units of ksTo. This also holds for the ground state analysis where 
the temperature is T = K but Tq = 298 K. Prom now on the pair elec- 
trostatic interaction will be written in reduced units so that Eq. (5) reads 
Ucoui = ZiZj/r. 

The macroion and the counterions are confined in a spherical impenetrable 
cell of radius R. The macroion is held fixed and is located at the center of 
the cell. The colloid volume fraction fm is defined as a /R^. To avoid image 
charge complications, the permittivity is supposed to be identical within 
the whole cell (including the macroion) as well as outside the cell. Typical 
simulation parameters are gathered in Table 1. 
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3 Ground state analysis 

In this section, we focus on counterion distribution exclusively governed by 
energy minimization, i. e. T — OK. In such a case correlations are maximal 
and all the counterions lie on the macroion surface. This situation has the 
advantage to enable accurate computation of energy variations in processes 
such as overcharging and also to provide a clear description of effects which 
are purely correlational in nature. The method employed here was successfully 
carried out in Refs. [12,13,16,17] and is explained in details in Ref. [17]. The 
Bjerrum length Is is set to lOo". Note that in the ground state the value of Z^, 
or equivalently the value of the dielectric constant e,., docs not influence at all 
the counterion structure. Only the electrostatic energy is rescaled accordingly. 



3.1 Neutral case 

First we consider the simple case where the system [macroion + counterions] 
is globally neutral. In order to characterize the two-dimensional counterion 
structure we compute the counterion correlation function (CCF) gdr) on the 
surface of the sphere defined as 



where c = Nd^na^ is the surface counterion concentration [Nc = Z^/Zc 
being the number of counterions) and r corresponds to the arc length on 
the sphere. Note that at zero temperature all equilibrium configurations are 
identical (except for degenerate ground state), thus only one is required to 
obtain gdr). The counterion correlation function gdr) is normalized as follows 



Because of the finite size and the topology of the sphere, g{r) has a cut-off at 
rgc — 7ra — 25.1a and g{rgc) = 0. Furthermore the absolute value of g{r) can 
not be directly compared to the one obtained with an infinite plane. 

Similarly, one can also define a surface macroion correlation function (MCF) 
gm{f^) for the microions (representing the colloidal structural charge) on the 
surface of the macroion. The normalization of gm{f) is very similar to Eq. (7) 




(6) 




(7) 
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and reads 

ira 

.j2.r,„Wdr=(Z„-l). 



(8) 



where the arc length has been rescaled by a factor a/ro so that gdr) and gm{i^) 
are directly comparable (see also the setup Fig. 1) and are defined in the same 
r range. 



3.1.1 Monovalent counterions 

We first treat the systems where we have monovalent counterions, that is we 
have to deal with the symmetric discrete system (-1:+1). The counterion corre- 
lation functions gdr) are computed for a central macroion charge [denoted by 
5f^^^^(r)] and for discrete macroion charge distribution [denoted by g^'^'^^r)]. 
Results for three structural charges = 60, 180 and 360 are given in Figs. 
2(a), (b) and (c) respectively. For the continuous case (central charge) the 
counterion structure consists of a pseudo-Wigner crystal (WC) as was already 
found in Refs. [12,13,16,17]. Also the higher the absolute number of counteri- 
ons Nc (i. e. the concentration c) the higher the order of counterion structure 
for the continuous case [compare Fig. 2(a) with Fig. 2(c)]. 

It turns out that in the case of discrete colloidal charges the counterion distri- 
bution is strongly dictated by the colloidal charge distribution and especially 
for low macroion surface charge density (Tm {Z^ = 60) [see Fig. 2(a)]. For 
Zm = 60, g^J-"^'-^\r) and gmir) are almost identical. This indicates that each 
counterion is exactly associated with one DCC site. The ground state struc- 
ture for Zm — 60 is depicted in Fig. 3(a) where one clearly observes this ionic 
pairing. 

This strong ionic pair association can be easily explained in terms of local cor- 
relations. Let us consider the picture sketched in Fig. 4 which holds for strong 
ionic pairing, where a given dipole A (ionic pair made up of a counterion and 
a DCC site) on the macroion surface essentially interacts with its first nearest 
surrounding dipoles B. Note that very similar lengths were also considered in 
a recent theoretical study [20] in the onc-dimcnsional case (counterion adsorp- 
tion on a linear polyelectrolyte) . It is important to have in mind that such a 
local description is physically justified due to the strong screening generated 
by ionic pairing. Thereby local correlations are twofold: (i) the attractive in- 
teraction between the DCC site of dipole A with its paired counterion and 
the counterions of dipoles B, and (ii) the repulsive interaction between the 
counterion of dipole A and counterions of dipoles B. The correlations between 
DCC sites are not relevant since they are fixed. The intra-dipole attractive 
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Fig. 2. Ground state surface correlation functions g{r) for monovalent counterions 
(Zc = 1) and for three macroion bare charges: (a) = 60 (b) Z^ = 180 and 
(c) Zm = 360. The two counterion correlation functions (CCF) gdr) are obtained 
for discrete colloidal charges [gcP'^'^\r) denoted by CCF^^^'^"-'] and for the central 
charge [gi^'~^\r) denoted by CCF^*^*-'^]. MCF stands for the discrete colloidal charges 
pair distribution gmif)- 

interaction Epm between the DCC site and its "pinned" counterion can be 
written as 



pm 



ZdZc 

a 



(9) 



For the elementary nearest inter-dipole (or inter-ionic pair) interactions, one 

can write for the attractive interaction between the DCC site of dipole 

A and the counterion of dipole B: 

= (10) 

CI'dc 
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Fig. 3. Ground state structures for discrete monovalent systems (-1:+1): (a) 
Zm = 60 and (b) = 360. The discrete colloidal charges (DCC) are in white, 
and the counterions in blue. Full ionic pairing association occurs. The correspond- 
ing counterion correlation functions gcij"^ can be found in Figs. 2(a) and (c). 

A similar expression can be written for the repulsive inter-dipole interaction 
involving counterions of dipole A and dipole B, which reads 

E^^ = ^. (11) 

(^cc 



Note that the repulsive counterion-counterion term E^^ alone, even if space 
truncated^, drives to the long-range ordered triangular WC structure. How- 
ever at zero temperature the DCC sites represent pinning centers for the coun- 
terions where the electrostatic potential is considerably lowered (due to the 

^ This statement holds if the cutoff is larger than the lattice constant. 



dipole A dipole B 




Fig. 4. Schematic view of the local electrostatic interactions and typical correlation 
lengths involved between nearest dipoles. The negatively charged DCC (-) are in 
grey and the positively charged counterions (+) in white. 
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term Epin) compared to its direct "empty" neighborhood (charge hole), which 
in turn prevents the counterions from adopting the ideal WC structure. This 
latter aspect was thoroughly discussed in Ref. [16]. Another important quan- 
tity characterizing discrete systems is the ratio 

Ppin ~1 (■^^) 
(^pin 



between the mean inter-dipole separation Sec (more exactly the mean counterion- 
counterion separation) and intra-dipole separation dpin of an ionic pair (in the 
present study dpm = a as depicted in Fig. 4). The value of Sec can be obtained 
by taking the first peak position of gi'"'^\r). 

Obviously, for sufficiently low macroion surface charge density am (i- e. large 
Ppin) the ionic pairing term Epin will be dominant and strong ionic pairing 
occurs. More specifically, when the typical inter-dipole distance is large com- 
pared to the intra-dipole distance then dipole-dipole interactions are weak (i. 
e., — E^^\ <^ \Epin\) and the DCC distribution dictates the counte- 

rion structure. This is what qualitatively explains our simulation findings for 
Zm = 60 [see Fig. 2(a) and Fig. 3(a)]. 

When (Tm becomes sufficiently important the situation may become quali- 
tatively different. In this case dipoles approach each other and because of 
excluded volumeQ Occ becomes comparable to a (see Fig. 4)[^. Thereby, the 
counterion-counterion repulsion term E^j^ (overcompensating -E+_) induces 
counterion ordering compatible with the local attractive pinning potential field 
generated by DCC centers. This effect can be inspected in Fig. 2(b) and Fig. 
2(c) where one sees that upon increasing am, g^'^'~^\r) is gradually less cor- 
related with gmir) and more correlated with g^'^\r). As a topological conse- 
quence, some counterions will be in contact with several (two or more) DCC 
attractors as can be seen in Fig. 3(b). 

The quasi-triangular counterion arrangement for high a^ (^m = 360) can be 
inspected in Fig. 3(b). For this symmetric system in size (same diameter for 
the counterions and the DCC ions) one expects that for a compact amorphous 
DCC layer the counterion structure should become perfectly ordered. This 
extreme limit which would correspond to unreachable experimental charge 
densities has not been addressed in our simulations. 

In parallel, increasing am induces by purely excluded volume effect a stronger 
local order within the DCC layer as can be checked on the gm{r) plots in Figs. 

^ Note that in the present model no surface dipole flip is allowed which should also 
be the case experimentally. 

^ The limiting case is where the global structure is compact, i. e. touching spherical 
microions. 
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2(a)-(c). This is quite similar to what occurs in a system of hard spheres where 
the (dense) hquid phase is locaUy correlated and the (dilute) gaseous phase is 
uncorrelated. 

In summary, the system depicted above is the siege of an order-disorder phase 
transition where upon increasing am (i- e. decreasing ppin) we pass from a 
disordered counterion structure (imposed by the DCC layer) to a long-range 
ordered one which is induced by local counterion-counterion correlations. 

Although results presented above concern one given random distribution (for 
each Zjn), we carefully checked that similar results and conclusions could be 
drawn for different random realizations (systematically five). This also holds 
for the following section below where we deal with multivalent counterions. 



3.1.2 Multivalent counterions 

We turn to the asymmetric discrete systems (—1 : +^c) where multivalent 
counterions are present {Zc > 1). The correlation functions g{r) for two 
macroion charges Z^, — 60 and Z^, — 180 and various counterion valences 
Zc can be found in Fig. 5. One remarks that upon decreasing the number 
of counterions Nc (i.e., increasing Zc) for fixed Z^, the first peak of g^r) is 
gradually shifted to the right (compare also the monovalent case given in Fig. 
2) whatever the nature of the macroion charge is (discrete or continuous). 
Furthermore, we observe for the discrete systems that upon increasing Z^ (for 
fixed Zm) the correlation between 

^(DCC)(^) 

and gm{i^) decreases and increases 
between gjP'~^'^\r) and gf^'^^r). This effect is clearly noticeable in Fig. 2(b) 
and Figs. 5(c-e) corresponding to Z^ = 180. The very high counterion valence 
Zc — 10 reported in Fig. 5(e) was undertaken in order to stress the counterion 
multivalence effect. These findings lead to the conclusion that the counterion 
valence has the effect of reducing the disorder in the counterion structure 
stemming from the randomness of the DCC distribution. 

This related phenomenon can be theoretically explained with simple ideas. 
Basically, the mechanisms involved in this counterion valence induced ordering 
stem from two concomitant sources: (i) topological and (ii) correlational. 

The topological aspect is due to the presence of {Z^ — Z^jZ^ unbound DCC 
sites (free of associated counterion) ensuring global electroneutrality [com- 
pare for instance Fig. 3(a) and Fig. 6]. It is to say that here, compared to 
the monovalent case (-1:-|-1), the counterions have all the more "freedom" to 
choose their pinning locations because Zc is high. To be more precise, the 
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Fig. 5. Ground state surface correlation functions for different multivalent systems: 
(a) Zm = 60, = 2; (b) = 60, = 3; (c) Z„ = 180, Z, = 2; (d) Zm = 180, 
Zc = 3; (e) Zm = 180, Zc = 10. The two counterion correlation functions (CCF) 
are obtained for discrete colloidal charges (DCC) and for the central charge (CC). 
The gm{f) curves (denoted by MCF) are identical from (a) to (b) and from (c) to 
(e). 
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Fig. 6. Ground state structure for (-1:+10) with Z^a = 180. The corresponding 
counterion correlation function gcij^ can be found in Figs. 5(e). 



number of topologically accessible "pinned" configurations is given by 





(13) 



m 



which reduces to 1 for = 1. In the ground state, counterions will "decide" 
to choose among these various possible arrangements the one which minimizes 
the total energy of the system. It is clear that this topological feature by itself 
promotes counterion valence induced ordering. 

Concomitantly, there is a purely counterion correlation induced ordering which 
is Zc dependent. Indeed, using similar arguments as those previously employed 
for monovalent systems (-1:+1) built on Eqs. (9-11), one can infer the role of 
Z^. More specifically, by assuming an ordered WC structure^ the term Ea^^ 
given by Eq. (11) can be rewritten as 



^ Rigorously, Eq. (13) holds when each counterion is associated with one and only 
one DCC site (case of low Om)- For high dm, it remains a good approximation to 
capture the essential physics. 

^ From a topological point of view, it consists in replacing the current (random) 
Voronoi structure by the ordered WC structure. 



E, 



'++ 




(14) 



a 
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where Occ in Eq. (11) is now given by 

a.. = See ~ ^ . (15) 



Equation (14) shows that for fixed Z„ and a (i. e., fixed macroion charge 

density) E^^^ ~ Z^/^ whereas i?pj„ ~ (recalhng that Z^ = I) and therefore 
for sufficiently high Z^. the term E^^' wiU be dominant. Thereby Z^. induces 
counterion ordering so as to minimize mutual countcrion-counterion rcpiilsion 
merely dictated by Eq. (14). As a topological consequence, some countcrions 
which would be in contact with several DCC sites if they were monovalent can 
now be in contact with less DCC sites (see Fig. 6). 

In summary, these discrete multivalent systems are again the siege of an order- 
disorder phase transition which is counterion valence controlled. 



3.2 Overcharge 



We now investigate the charge inversion (overcharge) phenomenon. The start- 
ing equilibrium configurations correspond to neutral ground states as were 
previously obtained. The method employed here is very similar to the one 
used in Refs. [12,16]. To produce a controlled overcharge, one adds succes- 
sively overcharging counterions (OC) in the vicinity of the macroion surface. 
Thereby the resulting system is no longer neutral. Using Wigner crystal con- 
cepts [6,21], we showed that the gain in electrostatic energy (compared to the 
neutral state) by overcharging a single uniformly charged macroion (i.e., cen- 
tral charge) with n overcharging counterions can be written in the following 
way [12,13,17]: 



AE^^ = A£;^°" + AE!^""" = ^ 



(7V, + n)3/2-7V3/2|+^2^. (16) 



2 



As before Nc — Z^/Zc is the number of counterions in the neutral state, A 
is the macroion area (47ra^) and « is a positive constant which was deter- 
mined by using simulation data for AE^^ . AE'^''\ which is equal to the first 
term of the right member, denotes the gain in energy due to ionic correla- 
tions. The derivation of this term can be found in Refs. [12,13,17], and the 
basic idea is that each counterion interacts essentially with its neutralizing 
uniformly charged Wigner-Seitz cell. The second term on the right hand side, 
AE^^"", is the self-energy of the excess of charge. This repulsive term stops the 
overcharging for sufficiently large n. Note that the WC concept for describing 
energy correlations is already excellent for highly short range ordered struc- 
tures (strongly correlated hquids, see Ref. [6] for a detailed discussion). The 
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total electrostatic energy of the system as a function of n is displayed Fig. 7 
(for monovalent counterions) and Fig. 8 (for multivalent counterions) for two 
bare charges = 60 and = 180. The energy curves corresponding to 
discrete systems were produced by systematically averaging over five random 
DCC reahzations. 

3.2.1 Monovalent counterions 

Let us first focus on the monovalent symmetric case (-1:+1) where for the 
neutral state each DCC site is exactly associated with one counterion as was 
shown above. The results in Figs. 7(a-b) show that the overcharging process 
occurring with a discrete macroion charge distribution is quite different from 
the one obtained with an uniform surface charge distribution. Especially for 
the smallest bare charge = 60, the effect of disorder is very important in 
agreement with what was already found above for the neutral state in Sec. 
3.1.1. The main effects of charge discretization are: (i) the reduction in gain 
of energy and (ii) the reduction of maximal (critical) number, n*, of stabi- 
lizing overcharging counterions (corresponding to a minimum in the energy 
curve). Both points were thoroughly discussed elsewhere [16] for an equivalent 
symmetric discrete system (-2:+2). It was shown that points (i) and (ii) can 
be explained in terms of ion-dipole interaction, which presently is the main 
driving force for overcharging. When the overcharging counterions are present, 
each of them will essentially interact (attractively) with its neighboring dipoles 
(ionic pairs). The attractive ion-dipole interaction increases with decreasing 
ion-dipole separation, i. e. increasing macroion charge density 0"^. This ex- 
plains why the energy gain increases with Z^ [compare Fig. 7(a) and Fig. 
7(b)]. On the other hand, the repulsion between the counterions is not fully 
minimized since they do not adopt the ideal WC structure that is obtained 
with a central charge which in turn explains (i) and (ii) . However for high bare 
charge {Z^ = 180) the overcharge curve obtained with DCC [see Fig. 7(b)] 
approaches the one from the continuous case as expected for high counterion 
concentration. This feature is fully consistent with what was already found 
in Sec. 3.1.1, where it was shown that the order of the counterion structure 
in the neutral state (for discrete systems) increases with am- In other terms, 
the WC approach through Eq. (16) is a good approximation for describing 
discrete systems at high since stronger ordering exists. 

Common features of overcharging between continuous and discrete systems are 
briefly given here. We note that n* increases with the macro ionic charge Zm- 
Furthermore, for a given n, the gain in energy always increases with Zm- Also, 
for a given macroionic charge Zm, the gain in energy between two successive 
overcharged states is decreasing with n. Note that at T = K, the value of 
acts only as a prefactor. All these features are captured by Eq. (16). 
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Fig. 7. Total electrostatic energy for monovalent counterions ground state configu- 
rations as a function of the number of overcharging counterions n: (a) Zm = 60 (b) 
Zm = 180. Overcharge curves were computed for discrete macroion charge distribu- 
tion (DCC) and macroion central charge (CC). The neutral case was chosen as the 
potential energy origin. Dashed lines were produced by using Eq. (16). For discrete 
systems (DCC) error bars are smaller than symbols. 

3.2.2 Multivalent counterions 



Now we are going to discuss the asymmetric discrete systems (—1 : +Zc) 
where multivalent counterions are present [Zc > 1). The results of figures 
8(a-d) indicate that the energy gain in the overcharging process at fixed Zm 
and n is higher the higher the counterion valence Z^ for both macroion charge 
distributions (discrete and continuous). For the continuous case this can be 
directly explained in terms of WC concepts [i. e. Eq. (16)]. Indeed the main 
leading term of the correlational energy AE!^°^ in Eq. (16) scales like 

AE^""^ ~ -Z^J^ (17) 



for fixed n and fixed macroion charge Zm, and recalling that Nc = Zm/Zc- 
Equation (17) quantitatively (qualitatively) explains why overcharging is stronger 
with increasing counterion valence Zc for the continuous (discrete) case. 

As far as discrete systems are concerned, the overcharging mechanisms occur- 
ring with multivalent counterions differ from those occurring with symmetric 
monovalent systems (-1:-|-1). This is again due to the presence of {Zm — Zm/Zc) 
unbound DCC sites in the neutral state as discussed in Sec. 3.1.2. When over- 
charging comes into play, each overcharging counterion becomes paired with 
some0 of these free DCC sites. Figure 8 shows that the overcharging with 

^ It can be one or more depending on the valence, surface charged density and the 
local DCC site arrangement. 
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Fig. 8. Total electrostatic energy for multivalent counterions ground state config- 
urations as a function of the number of overcharging counterions n: (a) = 60, 
Zc = 2 (b) Zm = 60, = 3 (c) Zm = 180, Z^ = 2 (d) Z„ = 180, = 3. Over- 
charge curves were computed for discrete macroion charge distribution (DCC) and 
macroion central charge (CC). The neutral case was chosen as the potential energy 
origin. Dashed lines were produced by using Eq. (16). For discrete systems (DCC) 
error bars are only indicated when larger than symbols. 



multivalent counterions (especially the higher Zm) is significantly less affected 
by colloidal charge discretization than in the monovalent case (see Fig. 7). 

For Zm = 60, simulations show that overcharging in the discrete case can 
even be stronger than in the continuous case [see Figs. 8(a) and (b)]. This 
phenomenon can be qualitatively understood by referring to the very low 
macroion surface charge density limit, where the correlation term /S.E^^ in 
Eq. (16) becomes negligible compared to the ionic pairing term Epi^ given by 
Eq. (9). In this limiting situation, the energy gain by overcharging is approx- 
imatively given by —nZaZc/dpin so that full overcharging occurs where each 
monovalent DCC site is paired with one multivalent counterion. 
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For = 180, the overcharging curves for discrete and continuous distribu- 
tions are almost identical [see Figs. 8(c-d)]. This is consistent with what we 
already found in Sec. 3.1.2 for the counterion structure in the neutral state, 
where we showed that g^'^'^\r) approaches g^^'^^r) with increasing Z^.. How- 
ever, the agreement between discrete and continuous cases is even better for 
overcharging than for counterion structure [see the corresponding gdr) given 
in Figs. 5(c) and (d)]. This is due to the fact that, as previously mentioned, the 
WC approach [Eq. (16)] quantifying the energy gain by overcharging is already 
excellent for highly short-ranged ordered systems. Generally speaking, all the 
ordering mechanisms related in Sec. 3.1 for neutral discrete systems hold for 
the overcharging features: all causes leading to ordering enhance overcharging. 



4 Finite temperature 

In this part, the system is globally neutral and is brought to room temperature 
Tq. We are interested in determining the counterion distribution as well as the 
counterion motion within the counterion layer. The cell radius R is fixed to 
4:0a so that the macroion volume fraction fm has the finite value 8 x 10~^. 

4-1 Strong Coulomb coupling 

The Bjerrum length Ib is set to lOcr as previously in the ground state study 
Sec. 3.1. In this section we consider two macroion bare charges Zm (60 and 
180) and three counterion valences Zc (1, 2 and 3). A typical parameter for 
describing the Coulomb coupling strength is the so-called plasma parameter F 
[22] defined asT — IbZ^/gcc- For our simulation parameters, F ranges from 2.6 
(for Zm = 60 and Z^ = 1) up to 23.1 (for Z^ = 180 and = 3). Under these 
conditions, systems are still highly energy dominated so that at equilibrium 
almost all (if not all depending on Z^ and Z^) counterions lie in the vicinity of 
the macroion surface (strong condensation). Therefore for the strong Coulomb 
coupling regime it is suitable to focus on the counterion surface properties. 
In the following sections we are going to study surface counterion distribution 
and diffusion. 

4.1.1 Counterion distribution 

Like in the ground state analysis, we characterize the counterion distribution 
via its surface correlation function. At non zero temperature, correlation func- 
tions are computed by averaging Zli^j S{r' —ri)6{r" —rj) over 1000 independent 
equilibrium configurations which are statistically uncorrelated. 



18 



The results for monovalent counterions are depicted in Fig. 9(a) and Fig. 9(b) 
for Zm = 60 and Zm = 180 respectively. For both charges Z^ the counterion 
distributions are weakly affected by charge discretization and 

^(DCC)(^) 

and 

^(cc)j^^-j are almost identical. A closer look on Fig. 9 reveals that the agreement 
between discrete and continuous distributions is even better for high macroion 
charge density {Zm = 180) as expected. In fact for monovalent systems the 
pinning term Epin has its lowest magnitude so that, for sufficiently high am, 
the fluctuating intra-dipole separation becomes comparable to the inter-dipole 
separation and discretization effects (i. e. pinning) are canceled. These pinning 
and unpinning aspects will be addressed in more details in the next section 
4.1.2. As expected, the counterion positional order for discrete and continuous 
cases is much weaker at room temperature than in the ground state case 
(compare Fig. 9 and Fig. 2). 

The results for multivalent counterions are depicted in Fig. 10. We now find 
that the counterion distributions are strongly affected by charge discretization, 
and especially the higher Zc. This is in contrast with what was found in the 
ground state analysis Sec. 3.1.2 where no counterion motion occurs. This effect 
is of course due to the pinning (inhibition of large counterion motion) which 
is proportional to Z^.. 

Note that all the statements above hold for the particular finite temperature 
Tq. However the effect of finite temperature discussed here should hold, at 
least qualitatively, for a large temperature range. For very low temperature 
one should recover all ground state properties mentioned in Sec. 3.1. 




Fig. 9. Surface correlation functions at room temperature Tq for monovalent coun- 
terions. The two counterion correlation functions (CCF) gdr) are obtained for dis- 
crete colloidal charges (DCC) and for the central charge (CC): (a) Zm = 60 (b) 
Zm = 180. MCF stands for gm{r)- 
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Fig. 10. Surface correlation functions at room temperature Tq for multivalent coun- 
terions. The two counterion correlation functions (CCF) gdr) are obtained for dis- 
crete colloidal charges (DCC) and for the central charge (CC): (a) Zm = 60, Zc = 2 
(b) Z„ = 60, Zc = 3 (c) Zm = 180, Z^ = 2 (d) Zm = 180, Z^ = 3. MCF stands for 
9m{r). 

4- 1-2 Surface diffusion 

This section is devoted to answer the following question: do the counterions 
only oscillate around the DCC sites or do they have also a large translational 
motion over the sphere? 

To study this problem we introduce the following observable: 
1 r 

Ax\t,to) = — — / dt[x{t ) -x{to)]\ (18) 

t tQ J 



which is referred as the mean square displacement (MSD), where x{to) repre- 
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sents the position of a given counterion at time t = to (at equilibrium) and 
x{t,to) is its position at later time t. All particles lying within a distance 9.2a 
from the macroion center are radially projected on the macroion surface of 
radius a — Sa to give x{t,to). The root mean square displacement (RMSD) 
Ax{t, to) is defined as 

Ax{t,to) ^ ^fA^^(t~to). (19) 



Like for the surface correlation function, the RMSD is measured on the spher- 
ical surface (arc length) and it has a natural upper limit na. For the case of 
free counterions (i. e. macroion central charge without pinning) the RMSD 
Ax free reads 



TT^ - 4 

Ax free = a\ — — ^ 13.7a. (20) 



This quantity Ax free will be useful to refer to the "unpinned" state. 

The results for discrete systems are sketched in Fig. 11 for — 60 and 
Zm = 180. Monovalent counterions are free to move over the macroion surface 
for both bare charges Zm considered here. Moreover, our simulation data show 
that the counterions gradually become pinned with increasing Zc- All these 
features are captured by the Zc dependency of the pinning term i?pj„. For 
multivalent counterions, the degree of pinning increases with decreasing Zm- 
This is due to the fact that the discrete charges get closer from each other 
by increasing Zm so that a counterion jump from site to site is energetically 
less demanding. For the continuous case, we have checked that for the same 
parameters counterions always have a large lateral motion and move all over 
the sphere. 



= 60 
= 180 




12 3 4 

z 

c 



Fig. 11. Root mean square displacement (RMSD) as a function of counterion valence 
Zc for Zm = 60 and Zm = 180. Errors are smaller than symbols. 
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4-2 Moderate Coulomb coupling 



In this last part, the Bjerrum length corresponds to that of water at room 
temperature {Ib = 2a = 7.14 A). For this moderate Coulomb coupling counte- 
rions occupy all the cell volume. Clearly, the probability of finding counterions 
plainly outside the macroion surface is no more negligible (in contrast with 
the strong Coulomb coupling). The target quantity is the fraction P(r) of 
counterions lying within a distance r from the macroion center and is defined 
as 

P{r)=N{r)/N^ (21) 



with 



N{r) = J 'iTcrjcv{ri)dri 



(22) 



where c„(r) is the profile of the volume counterion concentration and N{r) is 
the so-called integrated charge. 



The results for = 60 and = 12 are sketched in Fig. 12(a) and Fig. 
12(b) respectively. For the highest charge. Fig. 12(a) shows that discretization 
effects are canceled for any counterion valence. On the other hand, for the 
small charge density case. Fig. 12(b) shows that discretization effects become 
important for multivalent counterions. In the present situation, the Coulomb 
coupling is five times less important than in the strong coupling case studied 
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Fig. 12. Counterion fraction within a distance r from the macroion center for 
different counterion valence Zc- (a) Zm = 60 (b) Zm = 12. Data were obtained for 
discrete macroion charge distribution (DCC) and macroion central charge (CC). 
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in Sec. 4.1. Therefore pinning effects can only be noticeable for sufficiently low 
a^n (here — 12) and multivalent counterions. 



5 Conclusion 



We have performed MD simulations within the framework of the primitive 
model to study the coupled effects of macroion charge discretization and coun- 
terion valence. The macroion bare charge is carried by monovalent microions 
randomly distributed over the colloidal surface. Different correlational regimes 
were considered: (i) ground state and (ii) finite temperature. 

Concerning the ground state analysis, we were interested in the countcrion 
structure in the neutral state and the overcharging phenomenon. We demon- 
strated that the order in the surface counterion structure (disorder in counte- 
rion structure induced by the discrete random macroion charge distribution) is 
increased (decreased) by increasing macroion surface charge density a^, and/or 
counterion valence Z^. For monovalent counterions, we showed that the ratio 
between the intra-ion pair (made up of a discrete colloidal surface ion and 
a counterion) distance and the mean distance between ion pairs is a funda- 
mental quantity to describe counterion ordering. When overcharge comes into 
play similar effects occur. More precisely, for sufficiently high charge density 
a^n the overcharging with monovalent as well as multivalent counterions is 
quantitatively the same as the one obtained in the continuous case. For low 
a^n, the overcharging with multivalent counterions can even be stronger in the 
discrete case than in the continuous case counterions. In contrast, for mono- 
valent counterions overcharging is always weaker than in the continuous case 
but approaches the latter with increasing am- 

In the finite temperature case, strong and moderate Coulomb couplings were 
addressed. In the strong Coulomb coupling, we showed that counterion pinning 
is very weak for monovalent counterions but it increases with increasing Z^. 
and decreasing am- This involves an increasing disorder in the surface counte- 
rion structure with increasing Zc and decreasing am- In the moderate Coulomb 
coupling corresponding to an aqueous situation, the volume counterion distri- 
bution is only affected for low a^. and multivalent counterions. 

A future work will address the presence of added salts as well as the case of 
two interacting macroions. 
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